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o 

^H ABSTRACT 

o 

fSl Context. Mixing processes in stars driven by composition gradients as a result of the Rayleigh-Taylor instability are not anticipated. 

They are supported only by hydrodynamic studies of stellar convection. We find that such mixing occurs below the bottom edge 

^^ of convection zones in our multidimensional hydrodynamic shell convection models. It operates at interfaces created by off'-center 

^2 nuclear burning, where less dense gas with higher mean molecular weight is located above denser gas with a lower mean molecular 

^^ weight. 

Aims. We discuss the mixing under various conditions with hydrodynamic convection models based on stellar evolutionary calcula- 
CO tions of the core helium flash in a 1 .25 M© star, the core carbon flash in a 9.3 M© star, and of oxygen burning shell in a star with a mass 

^H of 23 Mo . 

Methods. We simulate the hydrodynamic behavior of shell convection during various phases of stellar evolution with the Eulerian 
hydrodynamics code HERAKLES in two and three spatial dimensions. Initial models for this purpose are obtained by state-of-the-art 
stellar evolutionary codes, namely GARSTEC, STAREVOL, and TYCHO for the core helium flash, core carbon flash, and oxygen 
^^ shell burning, respectively. Most of our analysis is performed for two-dimensional hydrodynamic models of shell convection during 

the core helium flash at its peak covering approximately 250 convective turnover timescales or 1 .4 days of stellar evolution. 
Results. The mixing manifests itself in the form of overdense and cold fingers enriched with matter of higher mean molecular weight, 
originating from density fluctuations at the lower boundary of the convective shell, and "shooting" down into the core. The fingers are 
Q neither produced by overshooting from the convection zone nor by thermohaline mixing. Instead, they result from the Rayleigh-Taylor 

^H instability at the lower convection zone boundary due to a negative mean molecular weight gradient (molecular weight decreasing in 

^ the direction of gravity). They do not appear when the mean molecular weight gradient is positive. We call this process Rayleigh- 

^ Taylor finger instability mixing or RTFI mixing for short. 
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^^ 1, Introduction ing conclusions (Deupree 1996; Dearborn et al. 2006; Herwig 
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et al. 2006; Meakin & Amett 2007; Mocak et al. 2009). These 



CN The hydrodynamic approach to model certain phases of stellar simulations also considered the effects of composition gradients, 

CNj evolution, like e.g., a nuclear core flash, by numerically solv- ^ut without investigating their influence in much detail, except 

lO ing the Euler or Navier-Stokes equations is essentially built f^^ ^^^ ^^^^ of Meakin & Amett (2007) highlighting the impor- 

O upon first principles in physics, and thus (almost) parameter- ^^^^ ^f turbulent entrainment (Fernando 1991) at the bound- 

O free. This approach is advantageous compared to (ID) stellar ^^.j^^ of convection zones. 
^—1 evolutionary calculations when modelling phenomena related to 

'^ flow instabilities which are inherently multidimensional in na- Composition gradients play an important role during the 

.!-J ture. One of these instabilities is the Rayleigh-Taylor (RT) in- evolution of stars of all masses. They can instigate double- 

><^ stability (Rayleigh 1883; Taylor 1949), caused by the buoyancy difl^usive instabilities (Kippenhahn et al. 1980; Turner 1985) 

H force which acts on fluid/gas elements (blobs) of matter that are '•^•' "salt-fingering" or the occurrence of vibrational overstability 

•^ over-dense or over-light with respect to the matter surrounding known as semi-convection (Grossman & Taam 1996, ;and refer- 

{[jgjjj ences therein) and can trigger gravitational settling (Stancliffe 

One consequence of the RT instability in stellar evolution & Glebbeek 2008) thereby influence the evolution of a star. 

is the occurrence of convection (Kippenhahn & Weigert 1990; Chemical gradients are also likely responsible that a star be- 

Weiss et al. 2004), whose role for stellar evolution is even nowa- ^on^es a red giant (Stanchfl^e et al. 2009), and they lead to deep 

days still mainly known from ID hydrostatic stellar evolution- fixing in envelopes of low-mass red giants (Eggleton et al. 

ary calculations. However, convection is a genuine multidimen- ^^^^)- 

sional phenomenon, /.^.,its properties can be reliably inferred Stellar layers with a negative mean molecular weight gra- 

only on the basis of multidimensional hydrodynamic Simula- dient (i.e., the molecular weight yu decreases in the direction of 

tions. To this end several hydrodynamic studies of shell convec- gravity) can develop a shell where thermohaline mixing occurs 

tion in the cores of stars (hence sometimes called core convec- (Kippenhahn et al. 1980), which is driven by a double-diffusive 

tion) have been performed in recent years studies, with interest- instability. Thermohaline mixing can be regarded as a diffusion 
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process where blobs of gas with higher molecular weight sink 
down, while other material rises and mixes with neighboring gas, 
eventually providing further rising and sinking elements. Such 
a situation is encountered, e.g., during the off-center ignition of 
helium in the helium core during the core helium flash, where a 
convective shell forms above the ignition point which is continu- 
ously enriched by nuclear ashes with a higher molecular weight 
(mainly ^^C). This creates a negative mean molecular weight 
gradient at the base of the convection zone, which should cause 
thermohaline mixing within a thermal timescale (Kippenhahn 
et al. 1980). Thermohaline mixing is supposed to manifest itself 
in the form of finger- like structures (Turner 1985), which are in- 
deed observed, e.g., in oceanography Williams (1974); Gargett 
& Schmitt (1982). On the other hand, results from experimen- 
tal chemo-hydrodynamics often reveal instabilities that have a 
finger-like structure as well, but which are not regarded as ther- 
mohaline mixing. Well known in this context (De Wit 2004, 
2008) is viscous fingering where a less viscous fluid is displacing 
a more viscous one in porous media (Saffamn-Taylor instability), 
and density fingering driven by a Rayleigh-Taylor instability in 
systems, where a denser solution resides above a lighter one. 

In the following we discuss a mixing process characterized 
by sinking dense blobs of gas that carry matter of higher mean 
molecular weight yu rapidly and presumably adiabatically down- 
wards (in the direction of gravity) into matter of lower//, leaving 
behind trails which resemble finger-like structures. This process 
that has not been anticipated by stellar evolutionary calculations 
is neither the result of thermohaline mixing, nor viscous fin- 
gering, nor density fingering. The mixing process operates be- 
low the base of convection zones where the mean molecular 
weight gradient is negative, i.e., where "lighter" matter resides 
above "denser" matter. Hence, it does not occur below convec- 
tion zones where the mean molecular weight gradient is posi- 
tive. As it is caused by a Rayleigh-Taylor instability, we call it 
Rayleigh-Taylor Finger Instability mixing, or RTFI mixing for 
short. Intrinsic properties of the finger-like structures inferred 
mainly from our hydrodynamic simulations of shell convection 
during the core helium flash rule out an explanation within the 
framework of thermohaline mixing as the finger-like structures 
are not hotter, but colder than the ambient medium. This fact also 
supports findings of Kippenhahn et al. (1980), who estimated 
based on calculations of Thomas (1967) for the non-central core 
helium flash of a 1.3 M© star that the timescale for thermoha- 
line mixing (the time it takes diffusion to double the geometri- 
cal thickness of the partially mixed region below the convection 
zone) is of the order 10^ yrs, which is much too slow to be seen 
in our hydrodynamic simulations. It is neither expected that the 
timescale will drastically change (decrease) with improved stel- 
lar models nor during the flash as the conditions are in general 
very similar (Sweigart & Gross 1978). 

Contrary to the thermohaline instability, a dynamic mixing 
process like RTFI mixing is unknown in experimental chemo- 
hydrodynamics. It is unexpected as matter is mixed from a 
less dense layer into a denser one. However, similarly looking 
buoyancy-driven instabilities were indeed found in the labora- 
tory at the interface between a solution of lighter hydrochlo- 
ric acid (HCl) residing above a heavier solution of sodium hy- 
droxide (NaOH) oriented vertically in the gravitional field by 
Zalts et al. (2008). The reaction of these two solutions at their 
common interface leads to the production of salt (NaCl), and 
to over-dense (with respect to the horizontal average) fingers of 
NaCl extending into the NaOH solution as shown in a snapshot 
from the experiment of De Wit (2008) or Almarcha et al. (2010) 
in Fig. 1. Such systems resemble principle similarities with the 




Fig. 1. Buoyancy-driven deformation of an interface between a 
less dense solution of HCl (yellow) residing above a denser so- 
lution of NaOH (blue) in an external gravitational field (pointing 
downwards in the figure). The field of view is 3 by 3 cm, and 
the snapshot is taken few minutes after contact (Almarcha et al. 
2010). 

Table 1. Some characteristic properties of our initial models: to- 
tal mass M, stellar population, metal content Z, mass Mm, outer 
(inner) radius R^, and nuclear energy production rate L^, respec- 
tively. 



Model 



M 
L 
O 



M 

[Mo] 



Pop. 



[Mo] 



[W cm] [lO^Lo] 



1.25 
9.3 

23. 



0.02 

0.02 

-'0.02 



0.45 
0.94 

2.67 



1.2(0.2) 
1.0(0.3) 
1.0(0.3) 



' 1 
0.01 
1000 



stellar shell convection models where lighter gas of higher mean 
molecular weight is residing above denser gas of lower mean 
molecular weight. Proposed explanations for this kind of labora- 
tory mixing range from amplified local density fluctuations due 
to salt production, a local density decrease due to a temperature 
rise up to differences in the mass diffusion coefficients of the re- 
acting species (Zalts et al. 2008). Recent work seems to favor 
the latter one (Almarcha et al. 2010). However, mass diffusion 
in stars is too slow to be responsible for the mixing observed in 
our simulations (see e.g., IMichaud et al. (2008)). 

The paper is organized as follows. In Sect. 2 we discuss the 
input models for our simulations. Our hydrodynamic code and 
the hydrodynamic shell convection models together with the lo- 
cal and temporal properties of the RTFI mixing are described 
in Sect. 3. A possible theoretical explanation of RTFI mixing is 
discussed in Sect. 4, while the breaking of dynamic stability in 
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Fig. 2. Temperature T (solid) and mean molecular weight// (dot- 
ted) as a function of radius for the initial core helium flash model 
(a), the core carbon flash model (b), and the oxygen shell burn- 
ing model (c), respectively . In each of these panels the region of 
interest below the base of the convection zone is highlighted by 
the shaded vertical strip. 



stars is analyzed in Sect, 5. Finally, a summary of our findings is 
given in Sect. 6. 

2. Initial Data 

We used three initial models for our simulations (Tab. 1, Fig. 2), 
which correspond to the helium core of a 1.25 M© star during 
the peak (maximum core luminosity) of the core helium flash 
computed with the GARSTEC code (Weiss & Schlattl 2000, 
2007), the carbon-oxygen (C-0) core of a 9.3 M© star at peak 
of the core carbon flash computed with the STAREVOL code 
(Siess 2006), and the core of a 23 M© star during oxygen shell 
burning (Meakin & Arnett 2007) computed with the TYCHO 
code (Young & Arnett 2005). The thermodynamic structure of 
all three initial models is qualitatively similar exhibiting of an 
ofl'-center temperature maximum due to nuclear burning in a par- 
tially electron degenerate stellar core. Directly outside this tem- 
perature maximum the inner boundary of a convective shell is 
located which extends towards larger radii. 

The core helium flash occurs in low-mass stars (0.7 M© < 
M < 2.2 M©) due to an ofl'-center ignition of helium by the 
triple-CK reaction after central hydrogen exhaustion in the semi- 
degenerate helium core. The core carbon flash ensues after cen- 
tral helium exhaustion and is characterized by ofl'-center car- 
bon ignition in a semi-degenerate carbon core of a rather mas- 
sive star (9M© < M < 12 M©), the dominant nuclear re- 
actions being i2c(i2c, a)^^NQ and i2c(i2c, p)^^Nsi followed 
by ^^0(a, 7)^^Ne. Oxygen shell burning is typical for mas- 
sive stars (M > 9M©) close to core collapse, when oxygen 
and carbon burn in a convective shell above a silicon- sulfur rich 
semi-degenerate core possessing an ofl'-center temperature max- 
imum (Arnett 1994; Bazan & Arnett 1998; Asida & Arnett 2000; 
Meakin & Arnett 2006, 2007). 

During the core helium flash the convective shell is en- 
riched mainly by ^^C which results in a negative mean molec- 
ular weight gradient below its base (Fig. 2 a) ^ . The situation is 
similar during the core carbon flash, where the nuclear ash from 
carbon burning increases the mean molecular weight inside its 
convective shell relative to that of the unburned inner carbon 
core (Fig. 2 b). Oxygen shell burning, which follows oxygen core 
burning, causes a positive composition gradient, /.^.,the mean 
molecular weight increases in the direction of gravity at the bot- 
tom of the convective shell (Fig. 2 c). 



3. Hydrodynamic Simulations 

The numerical simulations were performed with a modified ver- 
sion of the hydrodynamic code Herakles (Mocak et al. 2008). 
The code employs the PPM reconstruction scheme (Colella & 
Woodward 1984), a Riemann solver for real gases according to 
Colella & Glaz (1984), and the consistent multi-fluid advection 
scheme of Plewa & Miiller (1999). Self-gravity, thermal trans- 
port and nuclear burning are included in the code. Nuclear reac- 
tion networks are generated using a tool of Thielemann (private 
communication). In Table 2 we summarize some characteristic 
parameters of our 2D and 3D simulations based on the initial 
models given in Tab. 1 . 

The core helium flash simulations were performed on an 
equidistant spherical polar grid in 2D (r,0) and 3D (r,0,(f)). 
While the angular grid of the axisymmetric 2D simulations 
hefl.2d.2 and hefl.2d.3 covered the full angular range, i.e., 0° < 



^ The gradient grows for roughly 10 000 years since the onset of the 
core helium flash. 
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Fig. 3. Snapshots of the relative angular fluctuations of the helium mass fraction, AXC^He) = 100 x [X^Uq) - (XC^He))^] / (Xe^He))^ 
for the 2D model hefl.2d.l taken at ^ = 4638 s (a), and the 3D model hefl.3d in a meridional plane with azimuthal angle = -30° 
(b) and = +30° (c) at ^ = 6000 s, respectively. {)e denotes the angular average at a given radius. 



Table 2. Some properties of the 2D and 3D simulations: num- 
ber of grid points in r (A^^), 6 (Nq), and (A^^) direction, radial 
grid resolution Ar, angular grid resolution A^ and A0 in theta 
and /7/z/-direction, and maximum evolutionary time tmax of the 
simulation, respectively. 



model 


K 


N, 


N^ 


Ar 


Ae 


A0 


t 




# 


# 
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[lO^cm] 


n 


n 


[lO^s] 


hefl.2d.l 


180 


90 


- 


5.55 


1.33 


- 


30 


hefl.2d.2 


270 


180 
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3.70 
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30 


hefl.2d.3 


360 


240 
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2.77 


0.75 
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120 
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180 


90 


90 


5.55 


1.33 


1.33 
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360 


180 
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1.95 


0.5 
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10 


oxfl.2d 


400 


320 
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1.75 


0.28 
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6 < 180°, the 2D simulation hefl.2d.l was restricted to the an- 
gular region 30° < 6 < 150°. The 3D simulation hefl.3d was 
performed within an angular wedge given by 30° < 6 < 150° 
and -60° < cp < +60°, respectively. This allowed us to simulate 
the 3D model with a reasonable amount of computational time 
using a radial and angular resolution comparable to that of the 
2D model hefl.2d.l for several convective turnover timescales 
(Fig. 3). In all core helium flash simulations the radial grid 
ranged from r = 2 x 10^ cm to r = 1.2 x 10^ cm. Boundary 
conditions were reflective in all directions except for simulations 
hefl.2d.l and hefl.3d, where we utilized periodic boundary con- 
ditions in angular direction(s) to avoid a numerical bias in case 
of wide angular convective structures. Abundance changes due 
to nuclear burning were described by a reaction network con- 
sisting of the four a-nuclei "^He, ^^C, ^^O, and ^^Ne coupled by 
seven reactions. 

The core carbon flash simulation cafl.2d was performed on a 
2D equidistant spherical polar grid (r, 6) covering a 90° angular 
wedge centered at the equator (6 = 90°), and a radial grid rang- 



ing from r = 3 X 10^ cm to r = 1 x 10^ cm (Fig. 4). Boundary 
conditions were reflective in radial direction and periodic in an- 
gular direction. Abundance changes due to nuclear burning were 
described by a reaction network consisting of ^H, "^He, ^^C, ^"^N, 

16q 20]^^^ 22]^^^ 23]^^^ ^^^ 24]y[g ^ouplcd by 17 rcactions. 

Oxygen shell burning simulation oxfl.2d was performed on 
a 2D equidistant spherical polar grid covering a 90° wedge cen- 
tered at the equator, and a radial grid ranging from r = 2x 10^ cm 
to r = 1x10^ cm (Fig. 4). Boundary conditions were reflective 
in radial direction and periodic in angular direction. Abundance 
changes due to nuclear burning were described by a reaction 
network consisting of neutrons, ^H, ^He, ^^C, ^^O, ^^Ne, ^^Na, 

24Mg, 28Si, 31p, 32s^ 343^ 35q^ 36^^, 38^^, 39k, ^OCa, ^^Ca, ^4X1, 

46Ti, 48Cr, ^^Cr, ^^Fe, ^"^Fe, and ^^Ni coupled by 76 reactions. 

In the following subsections, we will briefly discuss one by 
one the main properties of our hydrodynamic shell convection 
models of the core helium flash, the core carbon flash, and the 
oxygen shell burning. This provides the basis for an assessment 
of RTFI mixing and its importance in a wider context. 

3.1. Shell convection 

3.1.1. Core helium flash 

The hydrodynamic properties of shell convection during the core 
helium flash are described in detail in Mocak et al. (2008, 2009). 
Convection starts early (t < 1000 s) and quickly extends over 
the whole convectively unstable region as predicted by mixing 
length theory (MLT). In axisymmetry, the convection is charac- 
terized by fast and large circular vortices, while in 3D the con- 
vective flow is less ordered showing slower and smaller turbu- 
lent features. The relative fluctuations of density p, temperature 
r, and carbon ^^C with helium ^He in the convection zone are of 
the order of 10""^, 10""^, 10"^ and 10"^, respectively (Fig. 3). The 
convective velocities in our 3D model match those predicted by 
MLT quite well (|v| < 1 x 10^ cm s"^), whereas the velocities in 
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Fig. 4. Snapshots of the relative angular fluctuations of the mass fraction of ^He during the core helium flash for model hefl.2d.3 
at ^ ~ 12000 s (a), of ^^C during the core carbon flash for model cafl.2d at ^ ~ 682 s (b), and of ^^O during oxygen shell burning 
for model oxfl.2d at ^ ~ 940 s (c), respectively. The fluctuations are defined by AX^^N) = 100 x [XC^N) - (XC^NM I {X{^N))e, 
where ^N g {"^He, ^^C, ^^O}, and {)e denotes the angular average at a given radius. 



our 2D models exceed those by up to a factor of four. Turbulent 
entrainment (Meakin & Arnett 2007) leads to a growth of the 
width of the convection zone on a dynamic timescale in both 
the 2D and 3D models due to an exchange between the poten- 
tial energy contained in the stratified layers at the boundaries of 
the convection zone and the kinetic energy of the turbulent flow 
inside the convection zone. 



3.1 .2. Core carbon flash 

An analysis of the hydrodynamic properties of shell convection 
during the core carbon flash, based on our 2D model cafl.2d, 
shows that the convective flow is dominated by small circular 
vortices. The maximum angular averaged velocities inside the 
convection zone |v| ~ 4 x 10^ cm s"^ exceed those predicted by 
MLT vmlt ~ 1.5 X 10^ cm s"^ by about an order of magnitude. 
This is an unusual result, as in general, velocities inferred from 
2D hydrodynamic simulations of shell convection typically ex- 
ceed MLT velocities by a factor of five at most. The relative an- 
gular fluctuations of density p and temperature T are of the order 
of 10"^ in the convection zone, while those of the carbon mass 
fraction X(}^C) are of the order of 10"^ (Fig. 4). Whether turbu- 
lent entrainment is operative could not be decided, as the sim- 
ulated timescales are too short. However, we intend to address 
this issue elsewhere. 



3.1.3. Oxygen shell burning 

The hydrodynamic properties of shell convection during oxy- 
gen shell burning are discussed in detail in Meakin & Arnett 
(2007). The findings of these authors are confirmed by our 2D 
model oxfl.2d very well. The angular averaged amplitudes of 
the convective velocities inferred from model oxfl.2d are |v| ~ 
1x10^ cm s"^ Turbulent entrainment has been detected in model 



oxfld.2d, too, but we have not further analyzed this phenomenon; 
for more details see Meakin & Arnett (2007). The relative angu- 
lar fluctuations of density p and temperature T are of the order 
of 10"^ in the convection zone, while for example those of the 
oxygen mass fraction X(^^O) are of the order of 10"^ (Fig. 4). 
The similarity of these values with those in our models of the 
core helium and core carbon flash implies that the fluctuations 
of thermodynamic quantities seem not to depend on the specific 
type of shell convection. 



3.2. Mixing below the convection zone In case of V^ < 

We discovered a new mixing process operating below shell con- 
vection zones during degenerate, off'-center helium and carbon 
burning (the core helium and carbon flash, respectively). 

The process is characterized by a steep negative mean molec- 
ular weight gradient V^ < with V^ = dln/j/dlnP, where 
lighter gas with a higher mean molecular weight resides above 
a denser gas of lower mean molecular weight. The (negative) 
gradient results from nuclear burning in the convection zone, 
and the assumption of instantaneous mixing used to calculate 
the corresponding stellar models. 

In the following, we present the thermodynamic character- 
istics of such a mixing process during the core helium flash in- 
ferred from our 2D model hefl.2d.3, which was simulated with 
the highest grid resolution. Mixing is caused by overdense blobs 
enriched by material from the convection zone (e.g., carbon ^^C) 
with higher mean molecular weight ji literally shooting down 
from the base of the convection zone towards the center of 
the star (/.^., in the direction of gravity). The blobs penetrate 
into denser layers of the core thereby creating finger-like struc- 
tures (Fig. 5 a, b) which we also find in our 3D core helium 
flash model hefl.3d (Fig. 5 c, d). The properties of the mixing 
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Fig. 5. Maps of the carbon mass fraction XC^C) (in units of 10"^) and its relative angular fluctuation AX(^^C) = 100 x (X(^^C) - 
{X(^^C)}o) I {X{^^C))e at the bottom of the convection zone in the 2D model hefl.2d.3 at ^ = 63430 s (two upper panels), and in a 
meridional plane of the 3D model hefl.3d at ^ = 6000 s (two lower panels), respectively. The vertical solid line marks the bottom 
boundary of the convection zone which is equal to the position of T^ax^ and the dashed line gives the location from where RTFI 
mixing is launched. {)e denotes the horizontal average at a given radius. 



found in our hydrodynamic core carbon flash model are similar. 
Therefore, we discuss them only briefly at the end of this section. 
A detailed analysis of the flow in the layers just beneath the 
base of the convection zone in model hefl.2d.3 revealed sink- 
ing dense blobs enriched by material with higher mean molecu- 
lar weight than that of the ambient matter into which they pen- 
etrate. Resulting finger-like structures become visible immedi- 
ately after convection starts in the layers below the temperature 
maximum. The finger-like structures can easily be spotted in 
maps showing relative angular fluctuations of the mean molecu- 
lar weight, density, and temperature, respectively (Fig. 6 a - c). ^ 
They grow continuously (Figs. 7 a), and reach a length of almost 



^ Note that in Figs. 5 and 6, which show the elongated finger-like 
structures, the maps are displayed in spherical coordinates (r and 6 in 
cm and degree) using a rectangular projection. 



2x10^ cm within the first 20 000 s of evolution. Their width re- 
mains almost constant at a value of roughly 1-2x10^ cm during 
the entire simulation. 

Although, the finger-like structures created by heavy sink- 
ing blobs of gas at earlier times eventually diminish ^, they are 
being constantly replaced by new heavy blobs again leaving be- 
hind trails of gas enriched by matter of larger yu. Thus, the lay- 
ers into which the fingers penetrate are enriched by material of 
higher//. Consequently, the initial value of yu = 4.06633 rises by 



^ The typical lifetime of a sinking heavy blob is about 350 s when our 
simulations have progressed for about 60 000 s. 



M. Mocak and E. Miiller: Rayleigh-Taylor finger instability mixing 




(e) 



135 



-0.4 



4.2 4.4 4.6 
r (10^ cm) 

Vr (10^ cm s" 




4.2 4.4 4.6 
r (10^ cm) 




(c) 



CD 





Ap/p (10-^) 




14U 




— '• 


120 


"^ 


100 


^^m 'M 


^ 


80 




1 


60 


'^^'m 


H 


40 


. . m 


^H 



0.8 
0.5 

0.3 

-0.0 

-0.3 

-0.5 
-0.8 



4.0 4.2 4.4 4.6 
r (10^ cm) 



^^^ Ledoux Criterion (10 ^) 



4.8 

5n 



CD 



140 
120 


^^^^^ 1 


n 


100 

80 

60 
40 


-^Hi 


f 




4.0 4.2 4.4 4.6 
r (10^ cm) 



4.8 



Fig. 6. Snapshots showing the base of the convection zone in the 2D model hefl.2d.3 at t = 63430 s. The upper four panels give 
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0.05% below the convective shell until the end of the simulation 
(Figs. 7a and 8) at a steadily increasing rate ^. 

The sink velocities of the overdense blobs are initially 
roughly 1 x 10^ cm s"^ reaching approximately 3 x 10^ cm s"^ 
at the end of the simulation at ^ = 33.3 hrs. The blobs are al- 
ways cooler and denser than the ambient matter by roughly 0.4% 
and 0.05%, respectively. Thus, the observed finger-like struc- 
tures cannot be caused by the double diff'usive instability which 
leads to "salt-fingers", as these are hotter than the surrounding 
matter leading to heat diff'usion that is fast compared to diff'usion 
of composition (Schmitt 1994). 

Using the size of relative density fluctuations Ap/p within 
the region where RTFI mixing takes place we can derive an an- 
alytic estimate for the velocity of the dense blobs causing the 
finger- like structures, if we assume that the blobs arise from a 
Rayleigh Taylor instability (Weiss et al. 2004). With Ap/p ~ 
5 X 10""^, inferred from our simulation (Fig. 6 a), we find 



'^A 



Ap 



7 X 10^ cm s" 



(1) 



where g ~ 10^ cm s"^ is the gravitational acceleration at the base 
of the convection zone, and A ~ 10^ cm the approximate length 
of the fingers (Fig. 7). 

The above estimate agrees well with the results of an anal- 
ysis of our simulation, where the gas inside the fingers sinks 
at a bit smaller velocities ranging from ~ 1 x 10^ cm s"^up 
to 3.5 X 10^ cm s"^ (Fig. 6 and 7). It confirms that the "buoy- 
ancy work" (~ v^) is consistent with the acceleration of gas 
seen in the simulation, and suggests that the mixing is the result 
of the Rayleigh-Taylor instability. Hence, we propose to call it 
Rayleigh-Taylor finger instability mixing (RTFI). The velocities 
observed inside the finger-like structures (Fig. 6 e) are smaller 
than those predicted by Eq. (1), likely because of the friction 
caused by our numerical scheme between the ambient medium 
and the sinking overdense blob. 

Note that the initial density fluctuations can result from wrin- 
kling of the contact discontinuity at the boundary of the con- 
vection zone due to the convective flow, turbulence, g-modes or 
p-modes. 

Eggleton et al. (2006) also found Rayleigh-Taylor insta- 
bilities driven by a molecular weight inversion caused by the 
^He(^He,2p)^He reaction in hydrodynamic simulations of the 
hydrogen-rich layers residing above the hydrogen burning shell 
in low-mass red giants. The corresponding mixing is qualita- 
tively similar to RTFI mixing. The mixing velocities observed 
by Eggleton et al. (2006) agree very well with those theoret- 
ically estimated for a RT instability according to the formula 
v^ ~ gHpA/d//u, where equal relative density and composi- 
tion fluctuations, i.e., Ap/p = Ayu/yu, were assumed. However, 
any density fluctuation will necessarily lead to composition and 
temperature fluctuations AT as well. Moreover, the formula of 
Eggleton et al. (2006) holds only for an ideal gas, and only if the 
blobs are in pressure equilibrium with the ambient medium and 
have the same temperature as the surrounding matter. 

The finger-like structures resulting from RTFI mixing in our 
simulations are in pressure equilibrium with the surrounding gas, 
as they move with subsonic velocities (the speed of sound is 
Cs ~ 10^ cm s"^ in the corresponding layers). However, they are 
always cooler than the ambient medium, i.e., T > Ti = T -\- AT, 
where T is the temperature of the ambient medium. Allowing 
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Fig. 8. Relative angular fluctuations of the mean molecular 
weight Ayu/yu = (jd-{ji}e) / {jd}e at the base of the convection zone 
in the 2D model hefl.2d.3 at h = 30000 s (a), t2 = 60000 s (b), 
and t3 = 120 000 s (c), respectively. In each panel the black solid 
line gives the angular averaged radial distribution of the mean 
molecular weight yu, and {)e denotes the angular average at a 
given radius. Vertical dotted line marks the initial position of the 
convection boundary at t = s as determined by Schwarzschild 
criterion. 
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Fig. 9. Snapshots from the 2D model hefl.2d.3 at t = 63430 s 
showing (a) the angular averaged density (solid) and the square 
of the Brunt- Vaisala frequency N^ (dotted), (b) the angular av- 
eraged temperature, and (c) the modulus of the relative den- 
sity fluctuation \(p - {p)e)l{p)e\ (including its angular variation; 
shaded) as a function of radius near the bottom of the convection 
zone. The vertical solid line marks the bottom of the convection 
zone (location of Tmax)^ while the dashed line gives the location 
from where RTFI mixing is launched. 



also for fluctuations of the mean molecular weight Ayu, pressure 
equilibrium (p = pi) implies for an ideal gas equation of state 



pT ^ {p + Ap)(T + Ar) 
jl jl -\- Afi 



(2) 



Solving for Ap/p this relation should provide an improved esti- 
mate for the relative density variation in Eq. (1): 



Ap 
P 



ji{T + AT) 



-1 = 



1 -h Ajilji 
1 -h AT IT 



-1. 



(3) 



However, in our models RTEI mixing operates in a partially 
electron-degenerate core, /.^.,the assumption of an ideal gas 
equation of state is not justified as it was the case in Eggleton 
et al. (2006). Using the formula of Eggleton et al. (2006) the ve- 
locity of the sinking overdense blobs of RTFI mixing can be es- 
timated to be V ~ 5 X 10^ cm s"^ {Ajilji ~ 3 x lO""^), while Eqs. 1 
and 3 lead to the considerably larger value v ~ 2x lO^cms"^ 
{Ajilji ~ 3 X 10-4, AT IT ~ -4 X 10"^), because the density con- 
trast inferred from Eq. 3 Ap/p = 4 x 10"^ is roughly an order of 
magnitude larger than that observed in the simulations. 

We point out here once more that although the finger-like 
structures resemble overshooting, they actually originate from 
density fluctuations just beneath the lower boundary of the con- 
vection zone (Fig. 9 a) where the gas is already cold compared 
to that in the convection zone (Fig. 9 b), but still partially en- 
riched by matter from the adjacent convection zone. The fingers 
are also neither directly connected with sinking plumes inside 
the convection zone nor with large density fluctuations seen at 
the lower convection zone boundary (Fig. 9 c). 

There exists a correlation ^ between the maximum value of 
the square of the Brunt- Vaisala frequency A^^^^ and the num- 
ber of new fingers (heavy blobs) appearing at a given time 



(Fig.7b). A^^ 



varies between 1 rad^ s ^ 



and 1.5 rad^ s ^ 



for 



the 2D model hefl.2d.3. This can be understood, as the fingers 
are caused by overdense blobs of gas and the Brunt- Vaisala fre- 
quency can be related to the amplitude of the density oscillations 
according to 



A^^ 



Hp \ 



ad ■ 



^„) 



Ap 
P 



N^ 



-Ar (4) 



Another well known interpretation of N^ relates it to the fre- 
quency of g-modes or internal gravity waves: 



A^^ 



Hp\ 



ad 



Ar = Aroe'' 



(5) 



where V^^^ = dlnTldlnP is the logarithmic temperature 
gradient in the ambient surrounding, Vad = (dlnTldlnP) ad 
gives the temperature change for an adiabatic displacement, 
and V^ = d In jl Id In P the corresponding logarithmic compo- 
sition gradient. Furthermore, we have 5 = -dlnp I dlnT and 
(p = dlnp I dlnji, where the partial derivatives are evaluated at 
constant values of (P,ji) and (P, T), respectively. Finally, Hp is 
the pressure scale height. 

The correlation tells us that larger density fluctuations (larger 
values of N^) result in the formation of denser and therefore 



^ The correlation is not very clear during the first 2 X 10^ s, as it is 
difficult to determine the number of fingers during this early epoch. To 
identify a finger one searches for overdense regions below the base of 
the convective shell that correspond to sinking blobs. Focusing on fluc- 
tuations in jl is misleading, as these often exist only as traces behind 
sinking heavy blobs. 
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heavier blobs of gas at the bottom edge of the convection zone, 
which sink more easily down into the core. Such blobs do not 
have to be extremely heavy (dense) initially. It is sufficient when 
their density remains higher than that of the ambient medium on 
length scales comparable with the size of the observed finger- 
like structures. The dense blobs sink down from the bottom edge 
of the convection zone to layers, which are about ~ 20% denser 
than those from where the blobs are created (Fig. 9). This is pos- 
sible as the sinking blobs are compressed, because they remain 
in pressure equilibrium with the ambient medium. Thus, the sit- 
uation is exactly opposite to that encountered in the case of con- 
vection, where hot blobs of relatively dense gas rise from the 
bottom of the convection zone over large distances because they 
expand, although their initial density is much higher than that of 
the layers where they eventually end up. 

No or only a weak correlation can be recognized between 
the flow in the convection zone and the occurrence of sinking 
dense blobs (causing the finger-like structures) when comparing 
the temporal evolution of the total kinetic energy of the convec- 
tive flow and the number of fingers appearing at a given time (up 
to six) or the maximum of the Brunt- Vaisala frequency squared 
(Fig. 7 c). However, there exists a correlation between the po- 
sitions of the hot spots at the base of the convection zone and 
the locations where the fingers appear. Fingers occur more fre- 
quently at the positions of the temperature maxima (by ~ 40%) 
as compared to any other place beneath the convection zone, the 
temperature maxima being typically connected with convective 
upflows. 

We have no explanation for the observed periodic behav- 
ior of N^ax^ ^^^ f^^ the occurrence period of the finger-like 
structures (Fig. 7 b), which is roughly 10^ s (or ~ 20 convective 
turnover timescales) at r ~ 60000 s and decreases at later times. 
The periodicity seems to be uncorrelated with the convection as 
the total kinetic energy of the convective flow and N^^^ are not 
clearly correlated (Fig. 7 c). 

We have also found RTFI mixing in multidimensional mod- 
els of shell convection during the core carbon flash, which sug- 
gested to us that a crucial ingredient for the operation of the 



instability is the presence of a composition gradient, not nec- 
essarily coinciding with the bottom of convection zone. During 
the core carbon flash mixing leads to the formation of several 
layers enriched by matter of a larger mean molecular weight 
separated by interfaces with high values of N^ (Fig. 10, a) indi- 
cating large density fluctuations (Eq. 4). The bottom two layers 
are completely detached from the convection zone whose lower 
edge is located at r ~ 4.6 x 10^ cm. The radial distribution of 
jd across these layers shows a staircase-like behavior, which is 
similar to that observed in experiments of thermohaline mixing, 
i.e., a thermohaline staircase (Krishnamurti 2003). 



3.3. Mixing below the convection zone in case of V^ > 

In order to test our preliminary conclusion based on 2D shell 
convection models of the core helium flash and the core car- 
bon flash, which both have a negative composition gradients 
at the base of the convection zone, we have simulated in addi- 
tion a 2D hydrodynamic model of shell convection during the 
oxygen shell burning in a massive star. The corresponding stel- 
lar model is characterized by a layer having a positive mean 
molecular weight gradient V^ > at the base of the convec- 
tion zone. For this situation we do not find the finger- like struc- 
tures of RTFI mixing beneath the base of the convection zone 
(4.15 X 10^ cm < r < 4.32 x 10^ cm) (Fig. 10b). 



4. Stability Versus Local Perturbations 

For assessing the dynamic stability of a stellar layer we consider 
a fluid element which is radially displaced by an amount Ar from 
its initial position (Kippenhahn & Weigert 1990). The density of 
the fluid element will difl'er at its new position from that of its 
surrounding by 



Dp 



ele 



Ar 



(6) 
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where {dpldr)eie is the density change of the fluid element when 
it rises a distance dr, and (dp/dr)sur is the density gradient of the 
surrounding matter. 

According to e.g., Kippenhahn & Weigert (1990) a stellar 
layer is unstable against convection, if a fluid element that is 
initially lighter than its surrounding will continue to rise when it 
is displaced against the direction of gravity (i.e., Ar > 0). Hence 
convection sets in when the following criterion is fulfilled: 



(a) 



^] -W <0, 



ad 



(7) 



where the differentials refer to changes in radial direction. If it is 
assumed that the fluid element moves sufficiently fast, it moves 
adiabatically, and hence (dp/dr)eie = (dp/dr)ad- Note that the 
same criterion also holds for a fluid element that is initially heav- 
ier than its surrounding, and hence sinks downward (i.e., Ar < 0) 
in the direction of gravity. Note further that such an over-dense 
fluid element will sink downward as long as its density remains 
higher than that of the surrounding matter. 

Allowing also for (radial) composition gradients within the 
stellar layers Eq. 7 leads to the well known Ledoux instability 
criterion (see e.g., , Kippenhahn & Weigert (1990)) 



^ad ■ 



JV,>0, 



(8) 



where V^^^, V^j, V^, and 0, S are defined above (see paragraph 
after Eq. 5). 

The Ledoux criterion determines dynamically unstable re- 
gions in a star. Hence, if we could identify layers in our hydrody- 
namic models beneath the base of the convection zone where the 
inequality Eq. 8 holds, we could prove that the over-dense blobs 
which start sinking downward from there and causing RTFI mix- 
ing, are the result of a Ledoux unstable stratification. Indeed we 



find that the left term of Eq. 8 i.e., V^ 



^ad 



rV^ is positive 



in our 2D core helium flash model in most grid zones where we 
observe RTFI mixing (see white regions in Fig. 6 f). 

As mentioned before RTFI mixing is similar to mixing by 
convection in stars characterized by hot blobs except that it is 
caused by sinking cold dense blobs enriched with gas of larger 
mean molecular weight. However, there exists another important 
difference between both dynamic mixing processes. In convec- 
tion the surrounding gas has to have a superadiabatic tempera- 
ture gradient in order for convection to set in, while RTFI mixing 
sets in even if the correspoding layers are initially dynamically 
stable and Ledoux instability criterion is not fullfilled. In the lat- 
ter case becomes the stability broken only subsequently in our 
simulations, when convection boundary becomes partially en- 
riched by material with higher jd from convection zone and its 
enriched dense blobs start to sink. 

As it seems that, the yu-gradient is the reason for the observed 
instability, we analyzed in more detail the dependence of dp/p 
on d/j/fi near the boundary of the convection zone. Using the 
equation of state p = p(P, T, yu), we rewrite the relative (radial) 
variation of the density as 



dp dP dT dii 
p P T II 



(9) 



where the partial derivative a = ^Inp/^lnP is taken at con- 
stant T and yu. Layers are stable as long as their density gradient 
(dplp)sur does not exceed the density change experienced by a 
displaced fluid element (dplp)eie (Eq. 7), which we assumed to 
be equal to the adiabatic gradient (dp I p) ad- 
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Fig. 11. Angular averaged relative density fluctuations for adi- 
abatic displacements (dashed line) of a fluid element in the 2D 
oxygen shell burning model oxfl.2d as a function of radius at 
t = 847 s (upper panel), and in the 2D core helium flash model 
hefl.2d.3 at ^ = 63430 s (lower panel). The red shaded area shows 
relative density fluctuations given by dp/p = adP/P - SdT/T -\- 
(pd/i/n (see Eq. 9), while the green shaded area gives the same 
quantity, however assuming a yu-gradient of the same absolute 
size but of opposite sign. The width of the shaded areas indi- 
cates the variation of the relative density fluctuations with polar 
angle ^ at a given radius. The vertical solid line marks the bot- 
tom of the convection zone (location of T^ax)^ while the dotted 
line gives the location from where RTFI mixing is launched. 



The radial distributions of (dp/p)sur and (dp/p)ad are shown 
in Fig. 11(a) for our 2D oxygen shell burning model oxfl.2d 
near the base of the convection zone . We analyzed the sta- 
bility of these layers by calculating the total differential dp/p 
(Eq. 9), which depends on the pressure, temperature, and com- 
position gradients (dP/P, dT /T and Jyu/yu, respectively). Given 
the positive yu-gradient of the model, we find that (dp/p)sur is 
always smaller than (dp/p)ad at the boundary of the convec- 
tion zone implying dynamic stability. Assuming for the com- 
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position gradient the same absolute value, but the opposite sign, 
we find that (dp/p)sur exceeds (dp/p)ad below the lower edge 
of the convection zone, i.e., the corresponding stellar layers at 
r ^ 4.2 X 10^ cm become unstable. 

The positive yu-gradient therefore reduces (dp/p)sur, making 
the layer more stable. The eff'ect of a negative //-gradient is ex- 
actly opposite leading to an increase of (dp/p)sur and a less stable 
situation compared to that with a positive yu-gradient. 

The destabilizing eff'ect at the base of our core helium flash 
model is rather small (~ 1% ) due to the smallness of the ji- 
gradient, but apparently strong enough to instigate RTFI mix- 
ing. This can be infered from Fig. 11(b), where (dp/p)sur ex- 
ceeds its adiabatic counterpart below the convection zone for 
4.35 X 10^ cm < r < 4.5 x 10^ cm. 

Hence, the driving element behind RTFI mixing is a negative 
/i-gradient. 

5. Discussion 

Despite the observed existence of dynamic mixing in stars 
driven by chemical gradients, its inclusion into stellar evolu- 
tionary calculations imply non existing observational evidence. 
Motivated by mixing driven by the ^He inversion found in hydro- 
dynamic simulations of hydrogen layers in low-mass red giants 
by Eggleton et al. (2006), Denissenkov & Pinsonneault (2008) 
discovered that the time which stars spend in the bump while 
climbing the red-giant branch is significantly prolonged. 

These finding suggests that layers with a negative gradient 
of the mean molecular weight should not become unstable on 
dynamical timescales. However, we observe just the opposite in 
our hydrodynamic simulations, where Wsur - ^ad - f V^ > is a 
very common situation, as demonstrated by Fig. 12. It depicts the 
stability properties of individual grid zones in the temperature 
(y sur - ^ad) v^ composition gradient ^V^ plane. 

Initially, all grid zones above and below the convection zone 
have a sub-adiabatic temperature gradient (V^^^ - Vad < 0) and 
thus are stable, while for grid zones in the convective shell the 
temperature gradient exceeds the adiabatic one by 10"^ render- 
ing them unstable (Fig. 12a). After the onset of convection the 
situation changes considerably (Fig. 12b). Whereas the temper- 
ature gradient remains super-adiabatic in many grid zones in- 
side the convection zone, some of these zones become sub- 
adiabatic and their composition changes. This can be understood 
by realizing that convection consists of rising hot plumes (super- 
adiabatic) of gas enriched by nuclear ash from nuclear burning 
and sinking cold plumes (sub-adiabatic) mixed with gas from 
layers further out. Chemical species are mixed quite uniformly 
throughout the evolved convection zone,/. ^., the corresponding 
grid zones tend to gather at V^ ~ 0. On the other hand, some 
zones from below the base of the convection zone become dy- 
namically unstable as their V^^^^ - V«j - ^V^ is positive. These 
are the zones which participate in RTFI mixing. 



6. Summary 

We have performed multidimensional hydrodynamic simula- 
tions of shell convection based on three diff'erent initial models. 
These correspond to helium core at the peak of the core helium 
flash in a low-mass star, carbon core during the core carbon flash 
and oxygen burning shell above silicon- sulfur rich core of a mas- 
sive star. 

We find a mixing process under the base of convection zones 
manifesting itself as over-dense blobs sinking in the direction of 
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Fig. 12. Distribution of zones from the 2D model hefl.2d.3 in a 
temperature (V^^^^ - V^j) v^- composition gradient (0/^)V^ plane, 
where each computational zone is represented by a symbol. 
Zones inside, above and below the convection zone are marked 
by red diamonds, blue asterices, and black crosses, respectively. 
The upper panel (a) gives the distribution of zones at ^ = 123 s, 
and the lower one (b) at ^ = 63430 s, respectively. The flow in 
the zones above (below) the diagonal line is Ledoux unstable 
(stable). 



gravity and originating in the convection boundary, which lead 
to appearance of overdense fingers enriched by gas of higher 
mean molecular weight. This mixing is likely a result of the 
Rayleigh-Taylor instability operating due to the presence of neg- 
ative mean molecular weight gradient (V^ < 0) . This is mainly 
confirmed by approximate match between theoretical velocity 
estimate for an element of matter with excess density over the 
average of its surrounding and observed velocities in our 2D sim- 
ulation with highest resolution. The theoretical velocities based 
on observed amplitude of density fluctuations in the fingers 
Ap/p ~ 5 X 10"^ give us a value of v ~ 7 x 10^ cm s"\ which is 
close to the observed velocities reaching values 1 x 10^ cm s"^ < 
V < 3.5 X 10^ cm s"^ The fingers are denser and colder than 
their surrounding, which rules out their possible explanation as 
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the thermohaline mixing. The observed fingers lead to a mixing 
of material with higher mean molecular weight from convection 
boundary partially enriched by gas from above-lying convection 
zone deeper to stellar cores and we call it Rayleigh-Taylor finger 
instability mixing or the RTFI mixing. The sinking over-dense 
blobs which cause the fingers are literally shooting from the con- 
vection boundary downwards. The number of such fingers ap- 
pearing at given time is correlated with maximum of square of 
Brunt- Vaisala frequency in the convective boundary, as it is pro- 
portional to the amplitude of density fluctuations. If the regarded 
amplitude of N^ is higher in the convection boundary, the fluctu- 
ations there are stronger (heavier) too. This in turn increases the 
likelihood of having heavy blobs of gas which can sink and cre- 
ate the fingers. The frequency at which the fingers appear is not 
clearly correlated with above-lying convection. However, pre- 
ferred spots from where the fingers start are correlated with tem- 
perature maxima just above the convection boundary, which one 
can associate with up-flows in the convection zone. 

We find similar mixing below the base of the convection zone 
in the carbon core during the core carbon flash, too. However, 
we did not observe it below the base of the convection zone in 
our model of shell convection in oxygen burning shell. The only 
qualitative diflTerence of the oxygen burning shell model com- 
pared to the others is a positive mean molecular weight gradient 
V^ (i.e., ji increasing in direction of gravity) at the bottom of its 
convection zone. Hence, we conclude that the condition V^ < 
must be the driving element of RTFI mixing. We could confirm 
this also analytically for the oxygen shell burning model by ex- 
amining the dependence of density fluctuations on the yu-gradient 
at the lower convection zone boundary. 

Implications of RTFI mixing for stellar evolution have still 
to be explored. While it almost certainly will not influence the 
core helium flash considerably, it could, for instance, have some 
important influence on the properties of the deflagration flame 
during the core carbon flash. 
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